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ABSTRACT 

The  use  of  an  optimal  prefilter  prior  to  sampling  can 
lead  to  reduced  mean  square  error  after  reconstruction. 

The  form  of  the  best  prefilter  or  the  associated  local  weight- 
ing function  depends  upon  the  reconstruction  method  to  be 
used.  We  display  and  discuss  these  weighting  functions  for 
the  most  common  spline  reconstruction  methods. 
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The  purpose  of  this  note  is  to  indicate  that  the  use 
of  the  standard  delta  function  filter  for  sampling  is  sub- 
optimal  when  standard  reconstruction  techniques  such  as 
nearest  neighbor,  linear, or  cubic  spline  interpolation  are 
used  for  signal  approximation.  The  errors  induced  in  re- 
construction using  these  schemes  have  led  to  various  suggest- 
ed modifications  of  the  reconstruction  methods.  Such  modi- 
fications are  often  impractical  because  of  their  computa- 
tional complexity.  This  note  points  out  that  mean  square 
errors  can  be  minimized  by  the  use  of  a specific  prefilter, 
i.e.,  by  modifying  the  sampling  technique. 

The  use  of  prefiltering  to  reduce  reconstruction  error 
is  not  new.  For  example,  it  is  known  that  the  ( L ^ - ) 
optimal  prefilter  for  delta  sampling  followed  by  (sinx)/x- 
reconstruction  is  a low  pass  filter  [1].  In  the  case  where 
the  original  signal  is  band-limited,  exact  reconstruction 
is  possible.  However,  in  precision  processing  of  digitized 
images,  for  example,  as  well  as  in  other  signal  processing 
applications,  (sinx)/x  interpolation  is  often  considered 
too  costly  when  interpolated  values  are  needed  for  geometric 
correction  and  display  magnification  [2].  The  usual  inter- 
polation schemes  include  nearest  neighbor  approximation, 
linear  interpolation,  and  cubic  spline  interpolation,  which 
we  may  regard  as  spline  interpolation  schemes  of  orders 
zero,  one,  and  three  respectively.  The  substitution  of 
spline  interpolation  for  sinx/x  reconstruction  calls  for  a 
new  prefilter;  despite  this  fact,  most  imaging  systems 
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attempt  to  sample  using  a pure  low  pass  filter  followed  by 
delta-function  sampling.  Indeed,  restoration  theories 
attempt  to  correct  for  inaccuracies  in  this  "ideal"  sampling 
filter,  even  in  the  case  of  schemes  which  represent  ideal 
reconstructions  by  cubic  splines  [3]. 

The  optimal  prefilter  for  a given  spline  reconstruc- 
tion scheme  is  the  Fourier  transform  of  the  optimal  weight- 
ing function.  We  will  determine  the  optimal  weighting 
functions  for  the  three  aforementioned  spline  schemes 
assuming  infinitely  many  equispaced  knots  (sampling  points) 
in  one  dimension.  The  problem  is  equivalent  to  that  of 
finding  the  best  \-2  spline  approximation  to  a given  noise- 
free  function  with  fixed  knots  [4].  By  using  the  Lp  norm, 
we  are  led  to  a linear  relationship  between  the  sample 
values  { y ^ > used  for  reconstruction  and  the  initial  signal 
f(x).  This  relationship  is  represented  by  the  weighting 
f uncti on  g^  , i . e . , 


00 

g i ( x ) f ( x ) d x 

- co 


(1) 


We  have  assumed  infinitely  many  sample  points  so  as  to 
assure  that  the  gi  will  all  be  translates  of  a single  opti- 
mal weighting  function  g,  which  is  the  desired  solution. 
Extensions  to  the  cases  of  finitely  many  points,  higher 
dimensions,  and  other  linear  reconstruction  techniques  are 
straightforward,  based  on  the  same  standard  numerical 
analysis  approach  used  below,  once  the  reconstruction  method 
is  properly  formulated. 


1 


2 


Denote  by  A the  transformation  which  takes  a sequence 


of  sample  values  { y ^ } 


into  the  interpolating  spline 


(or  reconstructed  function)  s(x)  with  equispaced  knots 

00 

located  at  points  of  the  sequence  { x ^ } . Suppose  that 

i = -«> 

A restricted  to  square  summable  sequences  &2  1s  a bounded 
linear  transformation  into  l_2 ( R ) . Then  for  a given  func- 
tion  f e L2.  the  sequence  y = (y^  which  minimizes 
| | A-  - f | | 2 can  be  shown  to  be  a solution  to  the  "normal 
equation" 

(A*A)+  = A*f  (2) 

where  A*  is  the  adjoint  transformation  from  (l_2(R))*  to 
(«.2)*,  i.e.,  from  L 2 ( R ) to  SL^.  Provided  A*A  is  nonsingular, 
the  optimal  sample  y^  is  a projection  of  the  linear 
operator  (A*A)-1A*  operating  on  f,  and  thus  the  optimal 
weighting  function  g^  is  the  l_2  dual  element  representing 
that  functional.  This  establishes  equation  (1). 

By  energy  conservation  considerations, 

.00 

! gi(x)dx  = 1,  and  so  in  fact  g..  e L ^ ( R ) . Thus  g.j(x)  is 

- 00 

the  optimal  weighting  function  for  general  Loo(R)  functions 
for  minimizing  the  uniformly  averaged  mean  square  error. 
Finally,  if  A is  covarient  with  respect  to  rigid  left  and 
right  shifts  of  the  samples  yi  -*•  yi  + k represented  by 
s(x)  -*•  s(x  + xk),  then  the  weighting  functions  g^  are  indepen- 
dent of  i modulo  interval  shifts,  thus  yielding  a unique 
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optimal  weighting  function  g(x). 

For  spline  reconstruction,  we  have 

. 

s(x)  = A * s l (x)  , (3) 

* j =,oo 
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where  the  4>  ^ ( x ) are  basis  spline  functions  chosen  to  suit 
the  particular  application.  If  we  regard  A as  an  infinite 
matrix  whose  kth  column  is  the  function  <t>k(x),  -°°  < k < °°, 
and  the  variable  x represents  the  row  number,  then  A*  = A^, 


entry  of  A*A  is  the  value 


4^  ( x ) 4> j ( x )dx  . 


desired  weighting  function  gi  is  then  simply  the  ith  row  of 
the  matrix  ( A*A ) ^ A*  , which  is  the  spline  generated  by  the 
values  {b-.}°°  , where  (b.,)  = (A*A)'^. 

1K  k = -oo  1J 

It  is  customary  to  use  B-spline  basis  functions  for 
the  4>.j(x),  which  are  shifts  of  the  B-spline  basis  functions 
at  zero  (Figure  1).  Since  the  B-splines  are  symmetric  with 
finite  support,  the  resulting  infinite  matrix  (A*A)  is 
banded  with  constant  bands.  For  zero  order  splines,  A*A 
is  simply  the  identity  matrix.  For  linear  splines,  A*A  is 
infinite  tridiagonal,  with  entry  1/6  in  the  codiagonals, 
and  2/3  in  the  diagonal.  For  cubic  splines,  A*A  has  seven 
nonzero  bands,  whose  values,  i nci dental ly,  correspond  to 
B7(x.j),  i = -3,  -2,  -1,  0,  1,  2,  3,  where  By  is  the  B-spline 
of  order  7 centered  at  Xg. 

In  each  case,  A*A  is  invertible,  so  that  we  can  solve 
for  the  optimal  samples  {y^}  used  for  reconstruction.  How- 
ever, for  the  cubic  spline  case,  a special  adjustment  is 
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needed.  Because  the  cubic  B-spline  has  support  three  in- 
tervals wide,  the  value  of  the  reconstructed  spline  at  the 
sample  points  is  given  by  the  formula 

s ( x i ) = (l/6)y._1  + (2/3)y.  + (l/6)y.+r 

Although  these  values  {y^  are  useful  for  constructing  s(x), 
since  the  cubic  B-spline  is  easily  computed,  our  usual 
notion  of  sampling  and  interpolation  requires  that  the 
sample  y^  satisfy  s(x^)  = y^  . Indeed,  y^  should  be  close 
to  f(x.j),  but  not  necessarily  equal  due  to  the  prefiltering. 
This  latter  condition  on  the  {y i ) is  satisfied  by  the  values 

y = C ( A * A ) ” 1 A*f , 


where 


(4) 


C 


1/6 


2/3  1/6 


for  A obtained  from  the  cubic  B-spline.  The  weighting 
functions  are  obtained  from  the  rows  of  C(A*A)_1A*.  The 
spline  s(x)  can  be  reconstructed  from  these  modified  samples 
using  the  unique  bounded  cubic  basis  splines  4'1-(x)  satisfy- 
ing H'.j(x)  = from  the  formula  s(x)  = 1°^  7 ^ ( x ) . For 

infinitely  many  equispaced  points,  the  4^  are  all  translates 
of  Yg,  shown  in  Figure  2.  Even  though  the  4^  have  infinite 
support,  they  decay  rapidly,  so  that  the  sum,  to  machine 
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precision,  is  finite.  In  fact,  we  could  use  the  'f^'s  to  de- 
fine A,  in  which  case  the  adjustment  using  the  matrix  C is 
unnecessary.  However,  using  a cubic  B-spline  to  define  A 
assures  that  A*A  will  be  easy  to  compute  and  will  be  banded. 

In  each  case,  the  appropriate  operators  are  coVariant 
with  respect  to  shifts,  so  that  we  obtain  unique  optimal 
weighting  functions,  shown  in  Figure  3.  As  expected,  for 
nearest  neighbor  approximation  the  optimal  samples  are 
obtained  from  local  unweighted  block  averages.  The  optimal 

weighting  functions  for  linear  and  cubic  spline  interpola- 

* 

tion  are  more  surprising,  despite  the  fact  that  an  optimal 
prefilter  is  a familiar  concept.  The  functions  are  not 
spline  approximations  to  (sinx)/x,  since  they  decay  exponen- 
tially. The  oscillatory  behavior  results  from  inverting 
the  positive  matrix  A*A,  and  is  governed  more  by  the  con- 
vergence of  the  partial  quotients  of  a continued  fraction 
than  by  the  function  (sinx)/x  [5],  Because  the  functions 


in  Figure  3 decay  so  rapidly,  truncated  versions  with  a 
support  radius  of  two  or  three  intervals  yield  very  nearly 
the  same  optimal  samples,  which  is  often  a distinct  advant- 
age over  (sinx)/x  weighting. 

In  any  case,  the  observation  that  the  weighting  func- 
tions behave  qualitatively  like  (sinx)/x  is  consistent 
with  the  earlier  result  that  (sinx)/x  reconstruction  calls 
for  low-pass  prefiltering  --  i.e.,  a (sinx)/x  weighting 
function.  Specifically,  if  we  use  basis  functions 
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to  define  A in  equation  3,  then  we  may  solve  (2)  to  obtain 

the  classical  weighting  function  for  Shannon  reconstruction. 

Since  (sine  * sinc)(x)  = nsinc(x),  a fact  which  follows 

o 

from  the  Fourier  equivalent  rect  (x)  = rect(x),  we  obtain 
(A*A)  = hi,  I the  identity  matrix.  So  the  rows  of  (A*A)~^A* 
are  all  translates  of  the  classical  low -pass  weighting 
function 


9 ( x ) 


s i n ( it/  h ) x 

( TT  / h ) X 


Error  analysis  shows  that  substantial  reduction  in 
expected  mean  square  error  can  be  achieved  by  using  the 
appropriate  prefilter  from  Figure  3,  as  opposed  to  unfilter- 
ed delta  sampling.  For  example,  if  the  autocorrelation 
function  of  the  initial  signal  decays  exponentially,  and 
the  sampling  is  very  fine,  a 35%  reduction  in  total  ex- 
pected mean  square  error  can  b^  obtained  by  using  optimal 
sampling  prior  to  linear  reconstruction,  as  opposed  to 
normal  delta  sampling  [5]. 

Implementation  of  these  results  involves  construction 
of  the  optimal  sampling  weighting  functions  for  local  inte- 
gration, or  their  Fourier  transforms  for  spectral  filtering. 
If  one  is  resampling  finely  sampled  data  to  achieve  com- 
pression, and  expects  to  use  spline  reconstruction,  then 
appropriate  weights  must  be  assigned  to  the  values  in  the 
neighborhood  of  the  new  sampling  point.  One  can  obtain 
these  weights  by  resolving  the  normal  equations  for  quanti- 
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continuous  case. 

To  construct  the  linear  weighting  function  of  Figure 
3b,  one  need  only  find  the  linear  spline  interpolating  the 
values  of  the  central  row  of  (A*A)”"*.  Thus  one  must  invert 
an  infinite  tridiagonal  matrix,  which  is  an  easy  procedure 
using  the  recursion  formula: 

1 

b°  Cl/(cQ/ Ci)2-4 


2c. 


2c- 


bk+l  -(co/cl )bk  ' bk-l 


* I 


• t 

K; 

* # 


where  Cq  is  the  entry  in  the  diagonal,  and  ci  is  the  off 
diagonal  entry.  In  this  case,  Cq  = 2/3,  and  c-|  = 1/6. 

The  middle  row  of  the  inverse  matrix  is  {...,b  ^ ,bp  ,b-j  , . . . } . 
Because  of  overflow  difficulties,  in  general  one  first 
solves  for  z k = b^/bk+1  when  calculating  the  b^s  [6].  In- 
cidentally, constructing  the  cubic  basis  function  of 
Figure  2 involves  the  same  matrix  inversion  to  obtain 
second  moments  at  the  knots. 

For  the  optimal  sampling  cubic  weighting  function 
(Figure  3c),  an  infinite  banded  matrix  with  seven  nonzero 
bands  must  be  inverted.  Since  the  matrix  is  diagonally 
dominant,  finite  versions  of  the  matrix  are  easily  inverted 
These,  as  it  turns  out,  converge  rapidly  to  the  desired 
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infinite  version  as  the  size  is  increased,  due  to  the  rapid 
decay  of  the  inverse  matrix  values.  A more  general  in- 
version technique  for  this  case  would  certainly  be  welcomed. 
In  both  cases, however , to  reconstruct  the  appropriate  spline 
weighting  function,  one  need  only  store  the  few  significant 
nonzero  values  ( bg  , b-j  , b^  , b^  , . . . ) corresponding  to 
9 ( x 0 ) , 9(xi ) » • • • » 9 ( x 3 ) > • • • • Since  g is  the  appropriate 
spline  interpolating  these  values,  it  can  be  obtained  easily 
from  spline  interpolation  routines  and  a short  table  of 
previously  computed  values  for  g(x^). 

To  su mm arize,  we  are  recommending  the  substitution  of 
one  of  the  weighting  functions  of  Figure  3 for  the  standard 
delta  function  sampling  method,  whenever  practical.  Our 
analysis  and  construction  of  these  weighting  functions 
assumed  the  case  of  infinitely  many  equispaced  samples, 
which  permitted  the  presentation  of  the  theory  of  optimal 
prefiltering  using  the  classical  Banach  spaces  and 
Since  the  resulting  weighting  functions  decay  rapidly,  they 
provide  very  nearly  optimal  samples  for  the  finite  case 
everywhere  except  within  two  or  three  sample  points  of  the 
endpoints.  In  addition,  their  effectively  narrow  width 
permits  easy  direct  implementation  of  the  local  weighted 


sum  for  sampling.  Extension  of  this  theory  to  sampling  in 
several  dimensions  is  straightforward,  especially  if  the  re- 
construction basis  functions  are  separable  in  each  variable. 
Finally,  the  reconstructed  function  may  be  considerably 
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